Three Dimensional Simulation of Jet Formation in Collapsing Condensates 
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We numerically study the behavior of collapsing and exploding condensates using the parameters 
of the experiments by E.A. Donley et al. [Nature, 412, 295, (2001)]. Our studies are based on a 
full three-dimensional numerical solution of the Gross-Pitaevskii equation (GPE) including three 
body loss. We determine the three body loss rate from the number of remnant condensate atoms 
and collapse times and obtain only one possible value which fits with the experimental results. We 
then study the formation of jet atoms by interrupting the collapse and find very good agreement 
with the experiment. Furthermore we investigate the sensitivity of the jets to the initial conditions. 
According to our analysis the dynamics of the burst atoms is not described by the GPE with three 
body loss incorporated. 

PACS numbers: 03.75.Fi, 42.50.-p, 42.50.Ct 
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I. INTRODUCTION 

Most of the single particle properties of Bose-Einstein 
condensates (BEC) in dilute weakly interacting gases are 
well described by the Gross-Pitaevskii equation (GPE) 
(for reviews see 0,0, y|). The GPE is well suited for in- 
vestigating static as well as dynamic properties of a BEC 
and also allows to investigate the stability of BECs with 
attractive interactions in a trapping potential. Inelastic 
processes which only become important for large parti- 
cle densities are usually accounted for by adding cubic 
and quintic damping terms to the GPE which are be- 
lieved to properly describe two- and three-particle loss, 
respectively. 

In early experiments on BEC with attractive inter- 
actions |j] the scattering length was fixed and limited 
the size of these condensates to a number N CI |5j which 
for typical experimental parameters was on the order 
of iV cr w 1000. In these experiments a collapse of the 
condensate was induced stochastically while growing the 
BEC out of a thermal cloud of atoms. In contrast, more 
recent experiments by Donley et al. [|| allowed to de- 
terministically induce collapses of the condensate which 
revealed a dramatic behavior for a 85 Rb BEC. In these 
experiments the sign of the scattering length is changed 
from positive (repulsive interaction) to negative (attrac- 
tive interaction) values by an external magnetic field. 
This sudden change in the sign of the interaction leads to 
a series of collapses and explosions, a dynamical behavior 
which provides a very good opportunity for testing the 
applicability of the GPE. 

In particular the nature of the atoms emitted in a burst 
during the collapse of the condensate is not very well 
understood at present Q • While some numerical studies 
d, 0| indicate that these atoms are produced coherently 
and can be described by the GPE, other numerical |l0| 
and analytical approaches 0, 0, 0, Q, 0] conclude 
that these atoms are not described by the GPE alone. 

The numerical studies arrive at different conclusions 



mostly due to the choice of the three particle loss rate 
near the Feshbach resonance. The work of Saito and 
Ueda @ predicts a series of collapses and explosions dur- 
ing the experiment and correctly reproduces the burst 
atoms by choosing a three particle loss rate that is smaller 
than the one used by Adhikari ^(| where one big collapse 
determines most of the behavior of the condensate and 
the burst atoms are not reproduced correctly. Varyin 
the three particle loss rate with the scattering length 
to match the burst energies allows to get good quantita- 
tive agreement with the experimental data. 

In the analytical approach by Duine and Stoof [l^ elas- 
tic binary collisions in the BEC are suggested to cause 
the bursts and in the dynamics of collective excita- 
tions driven by the collapsing condensate are investigated 
analytically and found to explain the overall features of 
collapsing and exploding condensates. The effect of a 
molecular state close to threshold near a Feshbach reso- 
nance on the interaction between the atoms is explored 
in [Til IT3 . Il4j . In fact, such novel interaction mechanisms 
arising close to Feshbach resonances which are not con- 
tained in the GPE could yield the relatively large burst 
energies seen in the experiment 12]. 

In this paper we study numerically the full three di- 
mensional GPE with attractive interactions including 
three particle loss. We first choose the three-body loss 
rate K® to match the observed remnant condensate par- 
ticle numbers and collapse times. We find only one single 
value for K® which fits the experimental data and subse- 
quently use this value in our numerics. Then we investi- 
gate the time evolution of the collapsing condensate. We 
find jets of atoms whose particle number and sensitivity 
to initial conditions is in quantitative agreement with the 
experiment by Donley et al. [fjj . The results of our com- 
parison can be seen in Fig. ^ where we plot the number 
of jet atoms as a function of the time t cvo i vc at which 
the collapse of the condensate is interrupted. Also, the 
shape and energies of the jets agree with the experiment. 
We do not, however, find burst atoms with energies as 
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FIG. 1: Number of atoms in the jets as a function of time 
^evolve [ms] (labelled as t) for aimt = 7ao, No = 16000 and 
^collapse = — 30ao. 'H — h': Experimental data from |J. '-o-': 
Numerical results for parameters equivalent to those used in 
the experiment. 



measured in the experiment and thus conclude that they 
are not described by the physics contained in the GPE 
alone. 

For our numerical studies we use the time-splitting 
sine-spectral method (TSSP) recently introduced by Bao 
et al. for solving the damped GPE H3 • This method 
is explicit, unconditionally stable, time transversal invari- 
ant, and of spectral accuracy in space and second order 
accuracy in time. It thus yields reliable results even in 
the case of having sharply peaked wave functions during 
the collapse. 

The paper is organized as follows: In Sec. [H] we will 
introduce the GPE including loss terms and present the 
numerical method we use to solve it. Then, in Sec. 11111 
we first adjust the three particle loss rate to match the 
experimental results for the number of remnant particles 
and the collapse time. This is followed by a detailed 
comparison of numerical and experimental results on the 
jet atoms and bursts of atoms. Finally in Sec. IIVI we 
draw some conclusions. 



II. MODEL AND NUMERICAL METHOD 



A. Model 

We consider a one component BEC with varying scat- 
tering length that is described by the GPE including a 
damping term that accounts for inelastic interactions 
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dt \ 2m 

for times t > with the initial condition 
V>(x,t = 0) = Vo(x). 



(1) 



(2) 



Here V( x > is the macroscopic wave function of the 
condensate, and x = (x,y,z) T is the spatial coordi- 
nate. We assume the trapping potential to be harmonic 
V(x) = m (uj 2 x 2 + uj 2 y 2 + lu 2 z 2 ) /2 with u x , oj v , lu z the 
trapping frequencies in a;, y, and z direction respectively, 
and m the mass of the atoms. 

The macroscopic wave function at time t — 0, i.e. the 
initial data ?/>0i is normalized 
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|Vo(x)| 2 d 3 x = 7V , (3) 



where A^o is the total number of condensate particles at 
time t = 0. The two-body interaction between the atoms 
is given by U{t) = 4Trh 2 a s (t)/m with a s (t) the s-wave 
scattering length. In the experiment the time depen- 
dence of a s (t) (cf. Fig. |3J) 6] is realized by changing the 
magnitude of an external magnetic field near a Feshbach 
resonance of the 85 Rb atoms. 

In the case of a positive scattering length a s the GPE 
with g = has a stable ground state solution V'gs(x) 
while if the sign of a s (t) is changed from positive to neg- 
ative the GPE becomes focusing and does not have a 
stable ground state solution anymore, i.e., the conden- 
sate may collapse. As in the experiment we will start 
our simulations at positive scattering lengths a s and for 
most of the calculations assume the BEC initially to be 
in its ground state 0, 0] corresponding to the initial 
condition Eq. J5| with "0o( x ) = V'gs(x). 

Loss from the condensate is described by the term 
g(p) = h ( K K 2 {t)p + K 3 (t)p 2 ) where K 2 (t) is due to two- 
body dipolar loss and K^{t) accounts for the three-body 
recombination inelastic processes. We assume the effects 
of K 2 (t) to be negligible and set Ks(t) equal to K° when 
a s (t) = a co n apsc , and otherwise [fjISllS- Furthermore 
we assume K° to be proportional to dcoiiapse El • Un- 
der this assumption we choose K® to match the number 
of remnant atoms in the condensate after the collapse 
and the collapse times observed in the experiments. We 
also compare its value with experimental data in |23| . 



In this section we will first introduce the GPE used 
for describing the dynamics of a BEC with varying inter- 
action strength. Then we briefly present the numerical 
method for solving this three dimensional nonlinear par- 
tial differential equation. 



B. Dimensionless GPE 

In our numerical computations, we simulate and 
present numerical results based on the following di- 
mensionless GPE. We introduce t = ui z t, x = x/Zq, 
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^(x,t) = l 3 / 2 j,{x,t )/VN~ , ^o(x) = l* ,z M*)/VN~o, 
with /o = \J^I niujz the size of the harmonic oscillator 
ground state, the dimcnsionless time and length units 
respectively. Substituting into Eq. (Q, multiplying both 

sides by l^ 2 /TiLu z yJ and removing all ~, we obtain the 
following dimcnsionless GPE 
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FIG. 2: Example of a time dependence of the scattering length 
a a during the experiment |6j . A condensate with ground state 
wave function Vo(x) is prepared at Qjs — flinit- Then the 
collapse is induced by going to a negative scattering length for 
a time r cvo i vc and finally the time evolution of the remaining 
particles is studied at large positive values of a s . 
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FIG. 3: Normalization, condensate widths and central density 
\ip(0, 0, 0, t)\ 2 as functions of time r cvo i V o [ms] (labelled as t) 
for Oi a it = 7a , a C oiiap SC = — 30ao and A^o = 16000. 



with the initial condition 



^(x,t = 0) = Vo(x); 



(4) 
(5) 



where j3{t) = Ana s (t)N /l , g{p) = <5i(t)/? c JV /> + 
8 2 (t)f3 2 N 2 p 2 , Si(t) and S 2 (t) equal to 6° and S% re- 
spectively when a s (t) = a co ii a pso, and otherwise, 
V(x) = \{j 2 x x 2 + rfy 2 + z 2 ) with -f x = uj x /u z and 

7 y = w s /w 2 , /3 C = 47r|a co i lapsc |/Zo, 5? = (^w z /3 c ), 

#2 = #3/ (^#0- For = Ocoiiapsc < we there- 

fore have j3(t) = —(3 c Nq. The macroscopic wave function 
ifj in the dimensionlcss GPE (0J is now normalized to 1 
at time t = 0, i.e., 



|<Kx,0)f d d x= / d J x=l. (6) 
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We assume that 5® and are independent of No and 

^collapse 

ana- 

The fraction of particles Nn (t) in a volume f2 is given 

by 



N a (t)= / |V(x,i)| 2 ei 3 x, 
Jo 



(7) 



and therefore the number of condensate particles in a vol- 
ume fl is Nq Na(t). Due to the loss term g(|?/>| 2 ) the total 
number of particles AWi(*) = N(t) = N N R s(t) is 
time dependent, decaying as 

NoUt) = -No f 9 m^t)\ 2 M(^t)\ 2 rf 3 x< 0. (8) 

We do not further care for the particles lost from the BEC 
by inelastic collisions and concentrate on the dynamics 
of the remaining condensate particles. However, parti- 
cles incoherently lost from the condensate might still be 
observed in an experiment. We also introduce the mean 
width of the condensate in directions x,y, z (cf. Fig. |3J) 
which is computed from the wave function as 

a 2 a (t)= f a 2 |V(x,£)| 2 d 3 x/N(t), a = x,y,z. (9) 

JR 3 



C. Numerical method 

We use the time splitting sine-spectral method (TSSP) 
which is described in detail in 0, for solving 
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the damped GPE. This method for numerically solving 
Eq. is based on a time-splitting, where at every time 
step one solves 



dt 



2 Vl 



(10) 



followed by evolving 



= (v(x)+mm 2 -i 9 -^P-)r- cl:d 



The linear Schrodinger equation without external poten- 
tial Eq. HI Up can be discretized in space by the sine- 
spectral method and then sovled in time exactly when 
homogeneous Dirichlet boundary conditions are applied 
46]. For each fixed x, the ordinary differential equation 
(ODE) Eq. 1)11(1 can be integrated exactly, too [1J|. In 
fact, the TSSP for the GPE combines the advantages of 
the spectral method which yields spectral order accuracy 
in space and can generate simple and explicit numerical 
formulae for partial differential equations (PDEs) with 
constant coefficients when a proper orthogonal basis is 
chosen, and the split-step method which can decouple 
nonlinear PDEs, e.g. the GPE, into linear PDEs with 
constant coefficients and a nonlinear ODE which can usu- 
ally be solved analytically. The merit of the TSSP for 
solving the GPE is that it is explicit, unconditionally sta- 
ble, time reversible and time transverse invariant when 
the GPE is, of spectral order accuracy in space and of 
second order accuracy in time, conserves the total num- 
ber of particles for the GPE Q without damping, i.e. for 
g = 0, and exactly reproduces the decay rate of the to- 
tal number of particles for the GPE (|TJ with a linear 
damping term, i.e. g = ao > 0. 

It is well known that the three-dimensional cubic non- 
linear Schrodinger equation (NLS) (i.e. GPE without 
external potential) with sufficiently large attractive two- 
body interaction leads to finite-time collapse of the spa- 
tial condensate density [23, El H2, i.e. \^\ 2 becomes 
a delta distribution for some finite value of time, and 
afterwards the solution ceases to exist. This collapse 
mechanism is halted by the three-body recombination 
term, which starts to kick in when the density becomes 
too large locally and then acts to instantaneously reduce 
the density in an explosion-like process [23, [^J . We re- 
mark that the mathematical analysis of the GPE with 
loss terms Eq. Q is not well developed yet, a systematic 
study has just begun right now. 



III. RESULTS 

Before presenting the results of our numerical studies 
we adjust the three particle loss rate by fitting the nu- 
merical results to the experimental data for the remnant 
number of condensate atoms and the time of the collapse. 
Then we study the formation of jet atoms and study their 
sensitivity to initial conditions. Finally we investigate 



bursts of atoms and present results of the simulations 
that do not agree with the experiment. 



Experimental parameters and three particle 
loss rate 



To determine the three particle loss rate we use the 
experimental results for a collapsing condensate of 85 Rb 
particles measured in [(| where a condensate of 85 Rb, 
with a particle mass m = 1.406 x 10 _25 kg, trapping fre- 
quencies uj x — Lu y = 2tt x 17.5/s, and lo z — 2ir x 6.8/s, 
i.e. a cylindrically symmetric geometry were used. In 
the experiment at time t — a BEC with ^V to tai(0) = 
Ao = Ao A'(O) = 16000 atoms and a scattering length 
of a s (0) = 7ao (where ao is the Bohr radius) is pre- 
pared. Afterwards, by changing the external magnetic 
field, the scattering length is linearly ramped down to 
&s = ^-collapse = — 30q,o (which corresponds to C = 
4.7717 x 10~ 3 ) in time t = 0.1ms and held there for a 
time T e voive as schematically shown in Fig. [21 This pro- 
cess leads to a strongly attractive and unstable conden- 
sate which undergoes a sequence of collapses and explo- 
sions before finally a remnant condensate with particle 
number A 7 JJ anneait is left. The remnant condensate is then 
measured at a positive scattering length (cf . Fig. [21 . 

We use the same time dependence for the scattering 
length a s in our numerics to find the quintic damping 
term 6®- Assuming cubic damping to be negligible we 
set Si = and fit the parameter of the dominant quin- 
tic damping term 6® according to the experimental re- 
sults for a co n a p S0 = — 30ao, i.e. the number of remnant 
atoms in the condensate after the collapse and the col- 
lapse times observed in the experiments. We simulate the 
experiment and adjust the three particle loss rate to nu- 
merically obtain a number of remnant atoms that agrees 
with the experiment. Fig. 0] shows the number of atoms 
in the condensate for different values of <5§. From this 
figure, we find a unique solution for the loss rate 6® given 
by (5° = 1.3 x 10 -3 which corresponds to a three parti- 
cle loss rate of K° = S^u z f3 2 = 6.756 x lO" 27 [cm 6 /s]. 
This value for three particle loss is in agreement with 
the measured values shown in Fig. 2c of Ref. [2j| and 
we use it for all of the following computations. We will 
check agreement of the numerics with the experimental 
results for the remnant condensate particle number and 
the collapse times. 



1. Remnant condensate particles 

In Fig. we show the comparison between the experi- 
mental and our numerical results for the number of rem- 
nant condensate particles Af remna nt as a function of the 
time T evo i ve ■ The results are in quantitative agreement as 
can be seen from Fig. [SJ We fit our results for JV rBn maat 
to a smooth function of the form N xemnaiIl t (^evolve ) = 

^remnant 4~ C^0 — ^remnant) ' ex P((^collapsc — r cvolvc) Adecay) 
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FIG. 4: Number of remaining atoms after collapsing a 85 Rb 
condensate of No — 16000 atoms with different values for the 
three-body loss rate 5°- Collapse is achieved by ramping the 
scattering length linearly from ai n i t = 7ao to a co n apsc = — 30ao 
in 0.1 [ms] as a function of time r evo i ve [ms] (labelled as t). The 
'*' are taken from the experiment 0- Curves are displayed in 
the order of decreasing iV° BIlmant for: Cyan: 5® = 0.00005 (i.e. 

= 2.598 x 10" 28 [cm 6 /s]); Green: 5§ = 0.00016 (i.e. k\ = 
8.315 x 10 -28 [cm 6 /s]); Blue: 8% = 0.0004 (i.e. K% = 2.079 x 
10" 27 [cm 6 /s]); Black: 5% = 0.0013 (i.e. K° = 6.756 x 10~ 27 
[cm 6 /s]); Red: <5§ = 0.003 (i.e. K$ = 1.559 x 10" 26 [cm 6 /s]); 
Yellow: 8° 2 = 0.008 (i.e. Kl = 4.157 x 10" 26 [cm 6 /s]). 



where iV r ° emnant gives the number of condensate parti- 
cles for Twelve — > do, tcoiiapse gives the time at which 
the condensate starts to collapse, and idccay is the de- 
cay time constant that determines the loss of particles 
during the collapse. The values found in the experiment 
were i C oiiapsc = 3.8, 8.6 [ms] and N° emnant = 5000, 7000 
for a G ollapse = — 30ao, —6.7ao, respectively which is in 
agreement with our numerical simulation where we get 
icoiiapsc = 3.6,9.35 [ms] and Ar° emnant = 5075,6970 for 
(icoiiapsc = — 30ao, — 6.7ao, respectively. We find decay 
times idccay = 2.8, 2.8 [ms] (cf. Fig. |3J| for a co i lapso = 
— 30ao, — 6.7ao, respectively, which also agrees with the 
experiment. Furthermore, from our numerical simula- 
tion for a C oiiapsc = -250a we find iV r ° cmnant w 1660, 
tdecay = 1.2 [ms] and icoiiapso = 1.1 [ms]. 



2. Collapse time 

Using the same value for S% we can also confirm the 
experimentally observed change in the time at which the 



condensate collapses as a function of the density of the 
initial condensate. In Fig. [S]we show the number of con- 
densate particles for Nq = 6000, a co ii a p SO = — 15ao and 
two different initial condensate densities with ai n ; t = 
and ai n it = 89ao. We find numerically i co iiapsc = 
5.7, 16.2 [ms] for a^it = 0, 89eto respectively, which is in 
excellent agreement with the experiment 6j. 



3. Nature of the collapse 

Our simulations reveal a series of collapses and explo- 
sions similarly to the experiment. In Fig. |3|we plot the 
particle density in the center of the trap as well as the 
widths of the condensate wave function in the different di- 
rections x, y, z. The fraction of condensate particle num- 
ber N(t) is also shown. The first collapse is marked by 
a sharp increase in the particle density in the center of 
the trap. During this first collapse a large amount of the 
particles is lost from the condensate. For the parameters 
chosen in Fig. [3] subsequent collapses are by far less im- 
portant than the first one as they have only a minor effect 
on the particle numbers. Also the peak density of these 
collapses is much smaller than for the first collapse and 
the widths of the condensate wave function are hardly 
affected. The times at which subsequent collapses hap- 
pen are determined by the stiffer oscillation frequencies 
in the harmonic trap (see Fig. ^J), i.e. as soon as those 
particles emitted during a collapse return to the z-axis 
the particle density in the center increases and the next 
collapse happens. 

The surface plots in Fig. [7] give a more detailed view 
of the evolution of the condensate density during the col- 
lapse. First the condensate contracts in the center of the 
trap and the particle density increases (Fig.[7^,b). Dur- 
ing the first collapse a number of sharp peaks forms in 
the vicinity of the trap center (Fig. 0J,d) which subse- 
quently, towards the end of the first collapse, spread out 
due to their kinetic energy (Fig. 0;,f). As we will show 
later these peaks are not of sufficient kinetic energy to 
produce the bursts of atoms seen in the experiment but 
interference effects between them lead to the forma- 
tion of the jets. 

Finally we note that our three dimensional simulations 
allow this degree of agreement with the experiment only 
for the three particle loss rate (5° chosen above. Simu- 
lations using smaller ^8] or larger [lCJ three particle loss 
rates do not yield numerical results that agree with the 
experimental data. 



B. Jet formation 

Next we are interested in the formation of jet atoms as 
observed in the experiment and thus simulate the follow- 
ing sequence for a s (t) (cf . Fig. |2J) : at t = the scattering 
length a s is ramped linearly from a- m it to a co n apse in 0.1 
[ms], then kept constant for a time r cvo i V c (during this 
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FIG. 5: Number of remaining atoms after collapsing a Rb 
condensate of iVo = 16000 atoms. Collapse is achieved by 
ramping the scattering length linearly from aj n i t = 7ao to 
acoiiapso = — 6.7ao,— 30ao and — 250oo in 0.1 [ms] as a func- 
tion of time involve [ms] (labelled as t). The '*' and 'o' are 
taken from the experiment H, the solid curves are our nu- 
merical solutions and the dashed curves are fitted to the ex- 
perimental points: N tots ,i(t) = N Tcrnnant (T cvo i vc ) 

^remnant 

(iVo ^remnant ) "6Xp( (^collapse devolve )/tdecay) with N° e 

7000, 5000, 1660; t co iiap S o = 8.6, 3.8, 1.1 [ms] and t dccay = 
2.8, 2.8, 1.2 [ms] for o co llapse = -6.7oo, -30a , -250a , 
respectively. 



time period we apply the damping term in our numerics). 
Then the collapse is interrupted by switching a s linearly 
back from <z C ollapse to a qucnc h = in 0.1 [ms], followed 
by changing a s exponentially from a que nch to a Gxpan d in 
5 [ms]. Then it is kept constant at a s — a oxpanc i. The 
resulting condensate density is shown in Fig. [S] where we 
can see the emergence of jet atoms and their dynamics. 

For Tcvoivc < ^collapse, i.e. before the condensate starts 
to collapse the outer region of the condensate is already 
affected by the negative scattering length and begins to 
expel some particles. As soon as the collapse has started 
this effect becomes more vigorous and condensate par- 
ticles are ejected from the core of the condensate dom- 
inantly in the xy plane forming jets. As the collapse 
continues, i.e. T OV oivc increases, the number of particles 
in the jet first becomes larger and then towards the end 
of the first collapse (cf. also Fig. [3] and Fig. for the 
duration of the first collapse) decreases again. Finally, 
when the collapse is allowed to complete no jets can 
be seen anymore. The numerical results for the num- 
ber of jet atoms as a function of r GVO i vc can be seen 
in Fig. ^ We count the number of atoms in the jets 
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FIG. 6: Number of remnant atoms in the Rb condensate for 
No — 6000 after ramping the scattering length linearly from 
abut = and ainit = 89ao at t = to a co ii ap sc = — 15ao at 
t = 0.1 [ms] as a function of time T eV oive [ms] (labelled as t). 



from the atom density IV'I 2 over a jet domain defined 
by [-5, -0.75] 2 x [-0.75, 0.75] U [0.75, 5] 2 x [-0.75, 0.75], 
and acquire the image 5.2 [ms] after applying a quenc h) 
|24|. These numerical results agree with the experimen- 
tal data very well and also the jet pictures in Fig.[H]which 
give an impression on the shape of the jets agree well with 
those observed experimentally. 

A large variance in the number of jet particles was 
found in the experiment. A possible reason for this could 
be that the initial condensate was not prepared exactly 
in the ground state. Then, if the jets are formed by in- 
terferences of particles emitted from different point like 
peaks in the atomic density Q as shown in Fig. [7] slight 
deviations in the initial condition from the ground state 
wave function may have a big influence on the numbers of 
particles in the jets. We have tested this assumption by 
introducing a small offsets of the initial condensate wave 
function from the center of the trap and indeed found 
large variations in the jet particle number. A typical ex- 
ample is given in Fig.[5]where an asymmetric jet of atoms 
can be seen. For small times t cvo i vc the behavior of the 
condensate is very similar to the case of a centered BEC 
shown in Fig.|S| This is in accordance with the small vari- 
ance in the particle number for r cvo i vc < 4ms observed 
experimentally. For larger times t cvo i vc the behavior of 
the condensate with an initial offset is significantly differ- 
ent from the centered case; an asymmetric jet is formed 
in Fig. 13 and we find large changes in the number of jet 
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FIG. 7: Surface plot of the density function \ip\ 2 in BEC with 
(Jinit = 7ao, ttcoiiapsc = — 30ao and No = 16000. At times 
Evolve (in [ms]) a). 0, b). 3.51, c). 10.53, d). 14.04, e). 
17.55, f). 21.06. 




FIG. 8: Jet images (i.e. image of \tp(x, 0, z, t)\ 2 ) for a series 
of t cvo i vc values for a in i t = 7a , iVo = 16000 and a co n apsc = 
— 30ao. The evolution times r cvo i vc were 2, 3, 4, 6, 8 and 10 
[ms] (from a) to f)). a quench = and a cxpand = 250a . The 
time from the application of a qucnc h until the acquisition of 
the images was equal to 5.2 [ms]. 



particles for small offsets in qualitative agreement with 
the experimentally observed values. The offset chosen in 
Fig. [5] corresponds to an initial potential energy of the 
condensate of E p w QA7Tuu z and thus to a temperature 
of T « 0.15nK. Therefore it seems plausible that even 
small temperatures - although they might not lead to an 
offset of the condensate as chosen in our numerical exam- 
ple - and/or uncertainties in the initial wave function can 
influence the jet formation and lead to the large variance 
seen in the experiment. 



C. Bursts of atoms 

Finally we also want to study the bursts of atoms 
observed in the experiments [(|. The bursts are parti- 
cles that are emitted from the condensate at relatively 
high energies when the collapse is allowed to complete. 
For finding the burst particles we compute <p z (z,t) and 
4>xy(x,y,t) as the axially and radially averaged density 



cross-sections, respectively, as follows: 

/oo />oo 
/ Wx,y,z,t)\ 2 dxdy, (12) 
-oo J — oo 
roc 

4>xy(x,y,t) = \ip{x,y,z,t)\ 2 dz. (13) 



Choosing a core domain Slo = [— b x ,b x ] x [—b y ,b y ] x 
[— b z ,b z ] and domains Q z = R \ [—b z ,b z ], fl xy = R 2 \ 
[— b x ,b x ] x [— by, by] [24[, we calculate the expectation 
value of the axial and radial potential energy per par- 
ticle in the region Q z and £l xy respectively 

Sa. 3*V*(*>*) dz 



^Axial(i) — 
^Radial (*) = 



Isk \ (J^ 2 + rfy 2 ) <Px V (x, y, t) dxdy 



fa xy <Pxy(x,V,t) dxdy 
and the number of atoms inside the core of the conden- 
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FIG. 9: Jet images for time r cvo ivc at 2, 3, 4, 6, 8 and 10 [ms] 
(from a) to f)) when we shift the center of the initial data 
from the origin to (0.375,0,0). All other parameters as in 

Fig.m 



sate 



N in (t) 



No I W 

/fin 



AUt(i) = N toUl (t)-N- m (t), 



where b x = b y = 1.5 and b z — 2.5 from our simulation for 
ainit = 0. Figure ITU1 shows these potential energies and 
the number of atoms in the condensate for No = 6000 
and a C oiiapse = — 30ao- We find qualitative agreement be- 
tween numerics and experiment in the number of burst 
atoms and the revival time, i.e. the time when the ra- 
dial energy i^Axiai(i) has a minimum. However, the en- 
ergy at which the atoms are emitted are far too small in 
our simulation and rather correspond to the energies at 
which the jet atoms are formed. Within our model of the 
three dimensional GPE with three particle loss we were 
only able to find particles emitted from the condensate 
at energies comparable to those observed experimentally 
if we decreased the three particle loss rate as done in 
the simulations of Decreasing the three particle loss 
rate leads to a more pronounced peak in the wave func- 
tion and higher particle densities during the collapse and 
thus to higher kinetic energies of the particles leaving 
the condensate. However, such small values of the three 
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FIG. 10: a) Radial and axial potential energy per particle 
outside the core of the condensate as functions of time r cvo i vc 
[ms] (labelled as r) and b) number of atoms inside and outside 
of the core of the condensate. We have chosen ai n it = 0, 
No — 6000, and <z C oiiapse = — 30ao. 



particle loss rate are ruled out by the above considera- 
tions in Sec. IIII Al Therefore we conclude that the three 
dimensional GPE with three particle loss is not able to 
describe all the features observed in the experiment 
correctly and some additional phy sical mechanisms like 
those investigated in pjl[l2l[l3 . ll5| need to be taken into 
account. 
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IV. CONCLUSIONS 

In conclusion we have shown that the GPE describes 
the physics of collapsing and exploding condensates apart 
from the energies of the burst atoms for the chosen three 
particle loss rate K°. We also found that no value for 
K® reproduces all the aspects of this experiment cor- 
rectly. We obtained excellent agreement for the number 
of remnant atoms and the collapse times. Also the jets 
of atoms are well reproduced by the GPE and we found 
that small variations in the initial condition for the wave 
function yield big changes in the number of jet atoms. 
The large fluctuations in the number of jet atoms ob- 
served in the experiment could thus be due to uncertain 
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